{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2a0b97d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import interp2d, interpn\n",
    "from numpy.matlib import repmat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00f6be88",
   "metadata": {},
   "source": [
    "# Entry/Exit with Homogeneous firms\n",
    "## and CES demand\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "30a07c92",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sequence_jacobian import solved, simple, create_model\n",
    "\n",
    "@simple\n",
    "def household(w, psi, nu):\n",
    "    L = (w/psi)**nu\n",
    "    return L\n",
    "\n",
    "@solved(unknowns={'v': 1.0, 'N': 1.0}, targets=['valfunc', 'accumulation'])\n",
    "def firms(C, sigma, Z, w, v, beta, delta, Ne, N):\n",
    "    mu      = sigma/(sigma-1)\n",
    "    p       = mu*w/Z\n",
    "    y       = p**(-sigma)*C\n",
    "    ell     = y/Z\n",
    "    pi      = p*y - w*ell\n",
    "    valfunc = v - (pi + beta*(1-delta)*v(1))\n",
    "    accumulation =  N - (1-delta)*N(-1) - Ne\n",
    "    return y, mu, p, ell, pi, valfunc, accumulation\n",
    "    \n",
    "@simple \n",
    "def mkt_clearing(L, N, Ne, ell, C, sigma, y, v, fE, w, Z):\n",
    "    labor_mkt = L - N*ell\n",
    "    goods_mkt = C- N**(sigma/(sigma-1))*y\n",
    "    free_entry = v - fE\n",
    "    return labor_mkt, goods_mkt, free_entry\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0a70a25",
   "metadata": {},
   "source": [
    "### Steady State"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d3fbc144",
   "metadata": {},
   "outputs": [],
   "source": [
    "ces = create_model([household, mkt_clearing, firms], name=\"CES Model\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bb885ca8",
   "metadata": {},
   "outputs": [],
   "source": [
    "calibration = {'beta':    0.96,\n",
    "               'psi':     1,\n",
    "               'fE':      1.0,\n",
    "                'nu':     2.0,\n",
    "                'delta':  0.11,\n",
    "                'sigma':  4.7924,\n",
    "                'w': 1,\n",
    "                'Z': 1.01, # note: things converge better if I set Z = 1.01. not sure why (or if it's ok to just do that...)\n",
    "                'Ne': 1.0, \n",
    "                'C': 1}\n",
    "\n",
    "ss = ces.steady_state(calibration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "283c6fca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.13520574056266665"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ss=ces.solve_steady_state(calibration, \n",
    "      unknowns={'w':1.0, 'Ne': .25}, \n",
    "      targets={'free_entry': 0.0, 'labor_mkt': 0.0})\n",
    "ss['Ne']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "21291b3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "ss=ces.solve_steady_state(calibration, \n",
    "      unknowns={'w':ss['w'], 'Ne': ss['Ne'], 'C': 6.5}, \n",
    "      targets={'free_entry': 0.0, 'labor_mkt': 0.0, 'goods_mkt': 0.0})\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Shock to fixed cost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4595522d",
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 300\n",
    "dZ = 0.5230094055387459*.01 * .685 ** np.arange(T)\n",
    "irf = {}\n",
    "irf = ces.solve_impulse_linear(ss, ['w', 'Ne', 'C'], ['free_entry', 'labor_mkt', 'goods_mkt'], {'fE': dZ})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0e8f721a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEYCAYAAABRMYxdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABHJUlEQVR4nO3dd3wc9Zn48c+zq94tS25ykY07tjFG2BSHZhJKCKRBICQBEkLIhVxy6UDuF+4ud8kluYQUEs4QCKHGBzgkhO7QQ7EN7gU3yZJc1Kzepef3x4zMWlZZ2VrN7uzzfr32pZ2yM89Iz46e+c53ZkRVMcYYY4zxk4DXARhjjDHGDDcrcIwxxhjjO1bgGGOMMcZ3rMAxxhhjjO9YgWOMMcYY37ECxxhjjDG+YwVOBIjIl0XkoIg0ishoj2IoFBEVkYR+pt8mIg+47ye7sQYHWea1IvJaJOI1ZigGy29jjImpAkdEikWkXUTyeo1f5+7sCj0KLTSWRODnwIdUNUNVq4dpuRErLlR1rxtrVySWHwnu33u613GYwcXC99YYr7gHlz2vbhFpCRm+2j0Y7eg133fcz74kIq3uuCoReVxExnu9TdEipgoc1x7gqp4BEZkPpHoXzlHGAinAZq8DMSaKDNv31lptjJ+4B5cZqpoB7AU+EjLuQXe2P4XOp6o/CVnETe5nZwI5wC9GdguiVywWOPcDnwsZvgb4Y+gMIvJhEXlXROpFpFREbguZliIiD4hItYjUishqERnrTrtWRHaLSIOI7BGRq/sKQESSReR2Ednnvm53x80Etruz1YrI3/v5/Gki8g93/etF5JyQaUfFICJzgDuB091KvXaw7QzxeTfG/SLyzX7iOaK5f7Dfg4j8TEQOudMuChn/koj80N22RhH5q4iMFpEH3RhXhx6ti8hsEXleRGpEZLuIXBEy7Q8icoeI/M2N4y0ROcGd9oo723p3PZ8SkTwRedL9ndaIyKsiEov57VcDfm8H+c725OcXRGQvcNT3SkQ+4bYUzXNz54ch084RkbKQ4WIRuVlEtrh5fK+IpAzz9hozolS1BngMmOd1LNEiFv8BvAlkicgccfqMfAp4oNc8TTg70xzgw8CXReSj7rRrgGxgEjAauBFoEZF04FfARaqaCZwBrOsnhluB04CFwEnAYuD7qvoecKI7T46qntf7gyJSAPwN+CGQC3wLeExE8vuLQVW3unG+4VbvOWFsZ49zgRnAh4Dvicj5/WxTT3yD/R6W4BRxecBPgN+LiIRMvxL4LFAAnAC8AdzrbutW4Ach63keeAgYg3N0/1sROTFkWVcB/waMAnYC/wmgqme5009yfx9/Ar4JlAH5OK1otwD2HJLoMdj3NpxcPhuYA1wQOlJErgP+GzhfVTeFGc/V7nJOwDny/f5QNsaYaCPOKeBPAO96HUu0iMUCB94/GvwgsA0oD52oqi+p6kZV7VbVDcDDODtHgA6cwma6qnap6lpVrXendQPzRCRVVferan+nma4G/l1VK1S1Euef8GfDjP0zwFOq+pQb3/PAGuDiIcYw2Hb2+DdVbVLVjTiFxlVHLehoA8VQoqp3uf117gPG4xQUPe5V1V2qWgc8DexS1RdUtRP4P+Bkd75LgGJVvVdVO1X1HZyjj0+GLOtxVX3b/eyDOAVlfzrcWKaoaoeqvqr2oLVo0+/3Nsxcvs3N5ZaQcV8Hvg2co6o7hxDLb1S11D3q/U/C+14Y45Ur3NbpnteEkGm/clv11wP7gW94EmEUiuUC59PAtfQ6PQUgIktE5EURqRSROpzWj7yQzz4LPOKeuvmJiCSqahPOUeWNwH731MjsftY/ASgJGS5xx4VjCnB5aLICS4HxQ4xhsO3sUTqUOMOI4UDIvM3u24yQ6QdD3rf0Mdwz7xRgSa/fw9XAuL7WBTT3Wk9vP8Vp5XnOPb32vQHmNd7o93t7DLnc49vAHapa1se0gQzpe2GMx1aoak7Ia1/ItH92xxWo6tXuQbchRgscVS3B6bR4MfB4H7M8BPwFmKSq2Tj9V8T9bIeq/puqzsU5/XIJbt8AVX1WVT+I0xKwDbirnxD24fyD7jHZHReOUuD+Xsmarqo/HiSGvloj+t3OEJOGGucQfg/HoxR4udfvIUNVv3wsC1PVBlX9pqpOAz4CfENElg1rxOa4DPK9DSeX+/oOfAj4voh8ImRcE5AWMjyOow35e2GMiS0xWeC4vgCc57Y49JYJ1Khqq4gsxjlqBEBEzhWR+W4/gHqcUxtdIjJWRC51+4a0AY1Af5dNP4yzU813z3v+P47uB9SfB4CPiMgFIhIUp9PzOSIycZAYDgITRSQpnO0M8a8ikub2bbkO+NNAwQ3x93A8ngRmishnRSTRfZ0qTofqcBwEpvUMiMglIjLd7Q9U78YcM5e9x5H+vrfh5HJfNgMXAneIyKXuuHXAxSKSKyLjcE5j9fYV9zuXi9Nfa8DvhTEm9sRsgeP281jTz+R/Av5dRBpwio8VIdPGAY/i/BPcCryMU3QEcDqq7gNqcM7//1M/y/8hTr+ZDcBG4B13XDhxlwKX4exUK3FaMr7trn+gGP6OszM/ICJVYWxnj5dxTt2sAn6mqs8NEuJQfg/HTFUbcI6+r3TXdQCno2hymIu4DbjPPb11BU5H6hdwCrI3gN+q6kvDHLY5TgN8b8PJ5f6WuR6nJfYuca7qux+nP0Ix8Bx9Fy8PudN2u6+wvr/GmNgh1g/TGBNPRKQYuF5VX/A6FmNM5MRsC44xxhhjTH8GLXBE5Ey3PwYi8hkR+bmITBnsc8YYY4wxXhn0FJWIbMC5md0CnHPbvwc+rqq971FhjDHGGBMVwjlF1eneMO0y4Jeq+kucKx6MMcYYY6JSOA+taxCRm3HuwHuWe3l1YmTD6lteXp4WFhZ6sWoTg9auXVulqvlex3E8LOfNUMV63lvOm6HqL+fDKXA+hXNPii+o6gERmYxz19gRV1hYyJo1/V0ZbsyRRKRk8Lmim+W8GapYz3vLeTNU/eV8OKeoGnBOTb0qztOyF+Lc6G44grpQnKdI77Rb65tYMFjOiuNX7vQNIrIoZFqxiGwUkXUiYntwExdsP2+8Ek4LzivAB0RkFM7N4tbgtOpcfTwrdk913YHz4L0yYLWI/EVVtxzPcv1CVelW56cCqqDuneoHu3VRvN3aSARSEoMjsJ6wcvYinJsOzsB58vrv3J89zlXVKoyJA7afH1hf+3lw9vW2nz9aatLQ9vPhFDiiqs0i8gXg16r6ExFZdyzB9bIY2KmquwFE5BGcjsz9Jn51dTXr1q1j4cKFdHV1cf/997No0SIWLFhAR0cHDz74IEVFRcybN4/W1lYeeeQRlixZwpw5c2hubmbFihWcfvrpzJo1i8bGRh599FGWLl3K9OnTqaurY+XKlZx11llMmzaNQ4cO8cQTT3DOOedQWFhIVVUVTz75JMuWLWPSpElUVFTw1FNPsej0s+lIyWZXSTm7332dwOSF1JFOc20lKQc3U54xizrSSGw5xKTWXWwKTqdWU8jpqmVW917e7p5GXXcyY6hlXnA/r7RPpUmTKAjUMT9hPy+3n0ALiUwK1HJiwgFebD+BNhKZEjjEnISDrGqfTgcJTA3WMCtYwfPtM+giyLRgNTODlTzbPhMlwPRgFdODVTzT7jw3c2awkqnBGp5tnwXA7GAFk4K1PN8+E4C5wYOMD9azqn0GAPMSDpAfaOTF9ukAzE/YT64083LHCQCclLCPbGnllQ7n6QknJ5STLu281jEVgFMSykiWTv7RUQjAqQmlBKWbNzucOw4sTtwLwNsdkwE4LbGELg2wutN5ZNAZicW0aQJrOycCsDRxD02axLudBcwcm8GNU6oZPXo0Z5/tXNz36KOPHnNiDiCcnL0M+KPbMf9NEckRkfGquj8SARkT5Xyxn1+w5AN0JOewu7Sc4nX/QCaeRJ1k0FJbSWrFZsoyZlHv7ucntu5is7ufz+48xGwt5a3uadSH7OdfbZ9Ko+3nh7SfT0kM8F8ntxy1nx83rq9HzTnCKnBE5HScFpsvuOOG43C5gCOf6FvGkUe6PSu/AbgBoKCgYBhWe+zaOrt5aXsFW1bXsntvOXm1Nfxm+2qqNJ1caWZxYhOrD5bSnZLNxKRWZnZ1o0BOWhKZycmkH0pg4fgcgmlZJLYECR48yCWFE0hMzUQbDtJVfojrZkwhkJJBV+0+2srq+NLsaQSTUumoKaOlvJ6b5k4nkJhCe9VeWsob+Pq8mQQSkmitKKZlXyPfnD8LCSbQenA3Lfub+M5Js5FAgJb9u2g50Mz3TnYSv3lfkLaKNr630B0uF9qqO/neAne4DNoOdXPqfGe4aW83HfXCknnucEknHY2HOP1EZ7ixuJ2u5nrOmOsO72mjq62JpbOd4YbdLWhHG2fNcod3NqHdXZwz0x3e0QjAeTOc4fr36pFAkGXT3eHttUhiMh+c5gzXbashmJzOBVNnk5uWBKXVEfqrHyGcnO1rngJgP87DIp8TEQX+V1WX97WS0JyfPHny8ERujDdiaj+vqrR0dPH8lgNsebOG3aX7GFdXw2+2rz1iP7+mohxSs5mQ2MLMrm4EGJWeRGZyChmHElg4IYeEtCwSmoMkHKzg0sIJJITs56+ZMYVgcjpddfvD2s9/bd4MggnJcb2fTwgIVK4b0t8znPvgnAV8C3hdVf9bRKYBX1fVfx7Smo5e7uXABap6vTv8WWCxqn61v88UFRXpSHc+K69t4Yl15fx1/X627q8HIDMlgQUTs5kxJpMTxmQwcVQqE7JTGZOZTHZqIoFA74cgGy+IyFpVLRrG5Q2asyLyN+BHqvqaO7wK+I6qrhWRCaq6T0TGAM8DX1XVVwZapxc5b2LbcOf98YiV/fzOikb+/G45f9u4nz1VznNgc9OT3P18BtPynf38+OwU8jNTyEpJwHmur4kG/eX8oC047g74lZDh3cBxFTeuMmBSyPBEnIcuRoW91c38+u87ePzdcrq6lUWTc/j+h+dw5vQ8Zo7NJGhFTDwKJ2f7nUdVe35WiMhKnOb7AQscY2JcVO/nN5XXcceLO3l60wECAqefMJovLJ3KmdPzKBydZkVMjAvnFFWkrAZmiMhUoBznqdKf9jAeADq6uvmf597j7ld3EwgInzt9Cp8/cyqTctO8Ds14L5yc/Qtwk9vXYAlQp6r73cedBFS1wX3/IeDfRzB2Y7wQlfv55vZO/uPJrTz89l4ykxP46nnT+dzpheRnJnsdmhlGnhU4qtopIjcBz+L06blHVTd7FQ9AaU0zX334XdaV1nJF0US++aFZjM1K8TIkE0X6y1kRudGdfifwFHAxsBNoBq5zPz4WWOkeESYAD6nqMyO8CcaMqGjcz28qr+OfH3mXPVVN3HDWNL5y7nSyUz25d62JMC9bcFDVp3D+IXhuXWktn/39WwD89upFXDx/vMcRmWjUV866hU3PewW+0sfnduM8082YuBJN+/nXdlTx+ftWk5uWxIPXL+GME/K8DslE0KAFjtu0+FWgMHR+Vb00cmGNrB0HG7j23rfJSUvkoetPs9NRxhjjM2/vqeGLf1zDtLx0HvriaeSmJ3kdkomwcFpw/ozzBPG/At0RjcYDpTXNfOb3b5EYDPDAF5ZYcWOMMT6zvrSWz/9hNeNzUnjg+iVW3MSJcAqcVlX9VcQj8UBHVzdf/OMaWtq7WHHj6UwZne51SMYYY4ZRU1snX3noHbJTnRb6vAzrSBwvwilwfikiPwCeA9p6RqrqOxGLaoQsf2U32w40cNfnipg9LsvrcIwxxgyznz67nbJDLaz40umMy7aLRuJJOAXOfOCzwHm8f4pK3eGYtbuykV+u2sGH54/ng3PHeh2OMcaYYbamuIb73ijmmtOnsHhqrtfhmBEWToHzMWCaqrZHOpiR0t2t3Pz4RlISAvzg0rleh2OMMWaYtXZ08Z3HNlCQk8p3LpztdTjGA4Ew5lkP5EQ4jhH11w37eGtPDbdcPIcxmdZkaYwxfvN/a0rZXdnEf31sPunJnt4RxXgknL/6WGCbiKzmyD44MXmZuKpy96t7OCE/nSuKJg3+AWOMMTGlu1v5/Wt7WDgphw/MsHvdxKtwCpwfRDyKEbS25BAby+v44Ufn2UMxjTHGh1Ztq6C4upnfXDDLnicVx8J52ObLIjIWONUd9baqVkQ2rMi55/U9ZKcm8vFFBV6HYowxJgLuenU3BTmpXHjiOK9DMR4atA+OiFwBvA1cDlwBvCUin4x0YJFQdqiZZzYd4KrFk0lLsnOyxhjjNxvKanl7Tw3XnVlIQjCcbqbGr8L5L38rcGpPq42I5AMvAI9GMrBIuP+NEkScJ4QbY4zxn9+/toeM5AQ+dar1sYx34ZS3gV6npKrD/FxUae/s5uG393LhieOYkJPqdTjGGGOGWXN7J89sOsDHFxWQmWJPCI934bTgPCMizwIPu8OfIkqeDDsUb++pob61k4+dbH1vjDHGj17bUUVbZzcXWN8bwyAFjjjdz3+F08F4KSDAclVdOQKxDatV2w6SnBDgzOl2yaAxxvjRqq0VZKYk2F2LDTBIgaOqKiJ/VtVTgMeHa6Ui8lPgI0A7sAu4TlVrh2v5vakqq7ZWcOb0PFKTgpFajTHGGI90dyurth3knFljSLTOxYbw+tK8KSKnDj7bkDwPzFPVBcB7wM3DvPwj7KpsZG9NM+fNHhPJ1RhjjPHIurJaqhrbOX+O7eeNI5wC51zgDRHZJSIbRGSjiGw4npWq6nOq2ukOvglMPJ7lDWbVVqeP9DJLfGOM8aUXthwkGBDOmWn7eePo9xSViExV1T3ARRGO4fPAnwaI4wbgBoDJkycf0wpWba1g7vgsxmfb1VPGGONHL2w9yOLCXLLT7Oop4xioBafnPjf3qGpJ79dgCxaRF0RkUx+vy0LmuRXoBB7sbzmqulxVi1S1KD8/P9ztOqy2uZ01JTXWemOMMT61t7qZ9w42cv7csV6HYqLIQJ2MAyLyA2CmiHyj90RV/flAC1bV8weaLiLXAJcAy1RVwwn2WLy0vZJuhWVzLPGNMcaPXtzudEOw/jcm1EAtOFcCrThFUGYfr2MmIhcC3wUuVdXm41nWYF7dUUVuehILCrIjuRpjjDEeeWfvIcZlpTBldLrXoZgo0m8LjqpuB/5bRDao6tPDvN7fAMnA8+6TXt9U1RuHeR0ArC+r5eRJOfbkcGOM8amNZXUsmGgHseZI4TxNfLiLG1R1+nAvsy+NbZ3sqmzkkgXjR2J1xhhjRlhdSwe7q5r4xCkRvRjXxCBf3w1pc3kdqnDSxByvQzHGGBMBm8rrAJhv3RBML74ucDaUuYlvTZfGGONL68tqAewUlTlKOA/bRETOAApD51fVP0YopmGzobyOgpxU8jKSvQ7FGGNMBGworWPK6DRy0pK8DsVEmUELHBG5HzgBWAd0uaMViPoCZ2NZrTVbGmOMj20oq+WUQnu4pjlaOC04RcDcSN6rJhLqmjsorm7m8qJJXodijDEmAiob2thX18rn7fSU6UM4fXA2AeMiHchw2+h2PLMOxsYY408by2sB62Bs+hZOC04esEVE3gbaekaq6qURi2oYbLDEN8YYX1tfWkdAYJ7t500fwilwbot0EJGwobSOwtFp9uA1Y4zxqQ1ltUwfk0F6cljXy5g4M+gpKlV9GdjG+49o2OqOi2oby+uYb6enjDHGl1SVDWV1LLD9vOnHoAWOiFwBvA1cDlwBvCUin4x0YMejqrGN8toWTrKOZ8YY40uVDW1UN7Uzb0KW16GYKBVOJ+NbgVNV9RpV/RywGPjXyIZ1fN470ADAnPGW+GZ4iciFIrJdRHaKyPf6mC4i8it3+gYRWRTuZ40x4Suudp7TPDU/w+NITLQKp8AJqGpFyHB1mJ/zTE/iF+bZk2XN8BGRIHAHcBEwF7hKROb2mu0iYIb7ugH43RA+a4xviMjlIrJZRLpFpGi4l19c3QTAlNy04V608YlwemY9IyLPAg+7w58CnopcSMevpLqJpIQA47NSvA7F+MtiYKeq7gYQkUeAy4AtIfNcBvzRvW/UmyKSIyLjce4EPthnj1JdXc26detYuHAhXV1d3H///ZQwhi1toxDtYnLNOxxKm0R96jgC3R1MOrSOmvTJNKSMJdjdzsRD66lOn0JjyhiCXW1MrN1AVcZUmpLzSOhqpaB2I1UZ02hKHk1iZzMT6jZTmXkCzUm5JHU2Mb5uCxWZM2hJyiG5o4Fx9ds4mDWT1sRskjvqGVe/nQNZs2hLzCKlo46x9e9xIGs2bYmZpLbXMqZhB/uz59KekE5aew35DbvYl30iHQlppLdVk9e4m/Kc+XQGU0hvqyKvcQ9lOQvoCiaT0VrB6KYSykadRFcgiczWg+Q27aV01EK6A4lktRxgVHMpe3MXoRIku2UfOc3llOSeAhIgu7mcnJZ9lIw+FYCc5jKyWg+wN9f5XzuqeS+ZrZXszT0FgNymEtLbaygddTIAo5uKSW2vpWzUQme4cQ8pnQ2U5ywAIK9xF0mdzezLmQ9AfsNOErtb2Zc9D4AxDTsIdrezP/tEAMbWb0fo5kDWHHd4GwAHs2YDMK5+K0qAg1mzABhft5muQBIVmTMAmFC3iY5ACpWZ05k7IYv5HdsZPXo0Z599NgCPPvoo48ZF1R09NgEfB/43EgvfW91MMCAUjEqNxOKND4TzNPFvi8gngDMBAZar6sqIR3YcSqqbmZybRiAgXodi/KUAKA0ZLgOWhDFPQZifBUBEbsBp/aGgoOD4IjbGI6q6FUAkMvvh4uomJo5KJTEY1ScUjIfCurZOVR8DHotwLMOmuLqJwtHWbGmGXV976t53+O5vnnA+64xUXQ4sBygqKtKFCxcCEAwGufbaa3vNvbTX8Fm9hs/uNXxOr+Fzew0v6zV8fq/hD/UavqDX8IW9hi8iep0e48MnHjH0yU9G9bUfAwot6idPnhzWZ3oOZI3pT7+lr4i85v5sEJH6kFeDiNQPx8pF5FsioiKSNxzLA+fSwZLqZqaMtv43ZtiVAaHP/pgI7AtznnA+a0xMEZEXRGRTH6/LhrIcVV2uqkWqWpSfnx/O/O6BrO3nTf/6bcFR1aXuz8xIrFhEJgEfBPYO53IrG9po6eiyFhwTCauBGSIyFSgHrgQ+3WuevwA3uX1slgB1qrpfRCrD+KwxMUVVezfxjYja5g4aWjuZYvt5M4Bw7oNzfzjjjsEvgO/QTzP9seq5gspacMxwU9VO4CbgWWArsEJVN4vIjSJyozvbU8BuYCdwF/BPA312hDfBGF/ouYLKWnDMQMLpg3PEiV4RSQBOOZ6VisilQLmqrh+sA9pQz81a4ptIUtWn6HUVoareGfJega+E+1lj/EpEPgb8GsgH/iYi61S1d6etY1Jy+EDWWnBM//otcETkZuAWIDWkz40A7bgdIAciIi/Q91PIb3WX27u3Yp96d7gcbP6S6iYSAsKEHLtE3BhjvOJebRuRK25LqpsRgUnWydgMYKA+OD8CfiQiP1LVm4e64P7OzYrIfGAq0NN6MxF4R0QWq+qBoa6nt+LqZiaOSiXBLh00xhhfKqluYnxWCimJQa9DMVEsnPvg3Cwio3DuzJoSMv6VY1mhqm4ExvQMi0gxUKSqVceyvN5Kqpus/40xxvhYcXUTk+30lBlEOJ2Mrwdewekc+W/uz9siG9axUVVKqprtCipjjPGxkupm62dpBhXOeZyvAacCJap6LnAyUDlcAahq4XC13hxq7qChrdNacIwxxqcaWjuobmq3/bwZVDgFTquqtgKISLKqbgNmRTasY3P4Cqo8a8Exxhg/siuoTLjCuUy8TERygD8Dz4vIIaL0DqwlPU+XtcreGGN8yQocE65wOhl/zH17m4i8CGQDz0Q0qmNUXNVMQGCiPV3WGGN8qaTGDmRNeAa6D06WqtaLSG7I6I3uzwygJqKRHYOS6iYm5KSSnGCXDhpjjB+VVDWTl5FERnJYz4o2cWygDHkIuARYy/tPQw79OS3i0Q1RSU2zNVsaY4yP7atroSDHWunN4Aa60d8l7s+pIxfO8TlQ18qZ04ftweTGGGOiTGVDGxNH2YGsGVw498F5QkSuEpGozqjubqWqsY38zGSvQzHGGBMhtp834QrnMvGfAx8AtorI/4nIJ0Uk6h70VNvSQUeXMsYS3xhjfKmzq5vqpnYrcExYwrmK6mXgZREJAucBXwTuAbIiHNuQVDa0AVjiG2OMT9U0taNq+3kTnrC6oYtIKvAR4FPAIuC+SAZ1LCoaWgEYkxl1jUvGGGOGQUXPgWyGFThmcIMWOCLyJ2AJzr1v7gBeUtXuSAc2VNaCY4wx/lbZ2LOfT/I4EhMLwmnBuRf4tKp2RTqY49FT2VsfHGOM8afDB7IZ1lJvBhdOJ+NXgJtFZDmAiMwQkUsiG9bQVdS3kZ4UJN1u/mSMMb7UU+DkWQuOCUM4Bc69QDtwhjtcBvwwYhEdo0q7dNAYY3ytsqGNjOQE0pLsQNYMLpwC5wRV/QnQAaCqLTh3M44qFfWt1sHYGGN8zA5kzVCEU+C0u1dRKYCInAC0He+KReSrIrJdRDaLyE+Od3mW+MYY42+VDW12BZUJWzjtfD/AuYJqkog8CJwJXHs8KxWRc4HLgAWq2iYiY45neQCV9W2cNcMS3xhj/KqqoY0546PqFmwmioVzo7/nReQd4DScU1NfU9Wq41zvl4Efq2qbu46K41lYS3sXDW2d1oJjjDE+VtnQxlkzbT9vwtPvKSoRWdTzAqYA+4F9wGR33PGYCXxARN4SkZdF5NQB4rhBRNaIyJrKyso+56m0S8SNMcbXWjvsQNYMzUAtOP/j/kwBioD1OC04C4C3gKUDLVhEXgDG9THpVne9o3BahU4FVojINFXV3jOr6nJgOUBRUdFR0yHkLsZZ1snYGGP8qNLuYmyGqN8CR1XPBRCRR4AbVHWjOzwP+NZgC1bV8/ubJiJfBh53C5q3RaQbyAP6bqIZhCW+Mcb42/t3Mbb9vAlPOFdRze4pbgBUdROw8DjX+2ecB3ciIjOBJOCY+/UcvotxliW+Mcb4kT2OxwxVOFdRbRWRu4EHcC4V/wyw9TjXew9wj4hswrmJ4DV9nZ4KV2VDG8GAkJtmd7c0xhg/sgLHDFU4Bc51OFc9fc0dfgX43fGsVFXbcQqlYVHR0EpeRhKBQNTdf9AYY8wwqGxoQwRy0+1A1oQnnMvEW4FfuK+oVNFgN/kzxhg/q2xsIzcticRgOD0rjAmvD07Uq2xos8c0GGOMj1XagawZIl8UOBUNbXYPHGOM8TErcMxQDVrgiMjl4YzzSle3Um3PoTLGGF+z51CZoQqnBefmMMd5orqpjW61uxgbY4xfqao9UNkMWb+djEXkIuBioEBEfhUyKQvojHRg4bJLB81IEJFc4E9AIVAMXKGqh/qY70Lgl0AQuFtVf+yOvw34Iu/fzPIWVX0q4oEb4wP1rZ20d3bbft4MyUAtOPuANUArsDbk9RfggsiHFp6KwwWOdTI2EfU9YJWqzgBWucNHEJEgcAdwETAXuEpE5obM8gtVXei+rLgxJkx2IGuOxUCPalgPrBeRlUCTqnbB4Z141GRZZb09aNOMiMuAc9z39wEvAd/tNc9iYKeq7obDjzm5DNgyMiEaEz1E5KfAR3Bu5roLuE5Va49lWfY4HnMswumD8xyQGjKcCrwQmXCGbnxOCh9eMN4qexNpY1V1P4D7c0wf8xQApSHDZe64HjeJyAYRuUdERvW3IhG5QUTWiMiayspjejybMdHgeWCeqi4A3uM4+m5mpiRw2cIJTMpNG7bgjP+FU+CkqGpjz4D7Pmqy7AMz8rnj04tISQx6HYqJcSLygohs6uN1WbiL6GNczyNIfgecgPMct/3A//S3EFVdrqpFqlqUn58/lE0wJmqo6nOq2tNf801g4rEua15BNr+88mQrcMyQhPOohiYRWaSq7wCIyClAS2TDMmbkqer5/U0TkYMiMl5V94vIeKCij9nKgEkhwxNx+rKhqgdDlnUX8OTwRG1MTPg8Tif9PonIDcANAJMnTx6pmIzPhVPgfB34PxHZ5w6PBz4VsYiMiU5/Aa4Bfuz+fKKPeVYDM0RkKlAOXAl8GqCnOHLn+xiwKeIRGxNhIvICMK6PSbeq6hPuPLfiXHn7YH/LUdXlwHKAoqKiY37wsjGhJJyHeItIIjALpwl+m6p2RDqwfuKoBEr6mJQHVI1wOJHkt+0Bb7ZpiqoOyzkeERkNrAAmA3uBy1W1RkQm4FwOfrE738XA7TiXid+jqv/pjr8f5/SU4lxm/qWQgmeg9cZLzoP/tsmr7Rm2vD9eInINcCOwTFWbw/yM5XzsiqqcH7TAEZEU4J+ApTg751eBO92HcEYFEVmjqkVexzFc/LY94M9t8pIff59+2ya/bc9QufeE+jlwtqoed295P/4+/bZN0bY94Zyi+iPQAPzaHb4KuB+Imsc1GGOMiTq/wbmlyPMiAvCmqt7obUgmnoRT4MxS1ZNChl8UkfWRCsgYY0zsU9XpXsdg4ls4l4m/KyKn9QyIyBLg9ciFdEyWex3AMPPb9oA/t8lLfvx9+m2b/LY9XvPj79Nv2xRV2xNOH5ytOB2M97qjJgNbgW5A3Zs4GWOMMcZEjXAKnCkDTVfVvnq7G2OMMcZ4JqzLxI0xxhhjYkk4fXCimohcKCLbRWSniBz1hOdoJyKTRORFEdkqIptF5Gvu+FwReV5Edrg/+312UTQSkaCIvCsiT7rDMb090cRyPjpZzkeO5Xx0ivacj+kCx32y+R3ARcBc4CoRmettVEPWCXxTVecApwFfcbfhe8AqVZ0BrHKHY8nXcPpq9Yj17YkKlvNRzXI+Aizno1pU53xMFzjAYmCnqu5W1XbgESDcByNGBVXd3/OcL1VtwEmWApztuM+d7T7go54EeAxEZCLwYeDukNExuz1RxnI+ClnOR5TlfBSKhZyP9QKnACgNGS5zx8UkESkETgbeAsb23Mrf/TnGw9CG6nbgOzhX2vWI5e2JJpbz0el2LOcjxXI+Ot1OlOd8rBc40se4mOw1LSIZwGPA11W13ut4jpWIXAJUqOpar2PxKcv5KGM5H3GW81EmVnI+nDsZR7MyYFLI8ERgXz/zRi33YaaPAQ+q6uPu6IM9T6AWkfFAhXcRDsmZwKXuQydTgCwReYDY3Z5oYzkffSznI8tyPvrERM7HegvOamCGiEwVkSTgSuAvHsc0JOI8pOX3wFZV/XnIpL8A17jvrwGeGOnYjoWq3qyqE1W1EOfv8XdV/Qwxuj1RyHI+yljOR5zlfJSJlZyP6RYcVe0UkZuAZ4EgcI+qbvY4rKE6E/gssFFE1rnjbgF+DKwQkS/g3EU61h9u6rft8YTlfEzx2/Z4wnI+pkTV9tiN/owxxhjjO7F+isoYY4wx5ihW4BhjjDHGd6zAMcYYY4zvWIFjjDHGGN+xAscYY4wxvmMFjjHGGGN8xwocY4wxxviOFTjGGGOM8R0rcIwxxhjjO1bgGGOMMcZ3rMAxxhhjjO9YgWOMMcYY37ECJ0qJiIrIdK/jMCaeiMgtInK313EYY46fFThDJCLFItIiIo0hr994HZfXROQ2EXnA6zjM8BORT4vIGjfX94vI0yKy1Ou4jpeInCMiZaHjVPW/VPV6r2Iy/iUiN4vIU73G7ehn3JUjG50/WYFzbD6iqhkhr5u8DsiYSBCRbwC3A/8FjAUmA78FLvMwLGNi0SvAmSISBBCRcUAisKjXuOnuvOY4WYEzTETkWhF5XUR+ISK1IrJbRM5wx5eKSIWIXBMy/x9E5E4ReV5EGkTkZRGZ0s+ys0XkjyJSKSIlIvJ9EQmISLKI1IjI/JB5x7gtTPk9R6gi8h13/ftF5KMicrGIvOd+9paQzwZE5HsisktEqkVkhYjkutMK3dNm14jIXhGpEpFb3WkXArcAn3KP8tdH6vdsRo6IZAP/DnxFVR9X1SZV7VDVv6rqt938u11E9rmv20Uk2f1sT+59MyT3rgtZ9sUissXN/XIR+ZY7/loRea1XHIdP17rfm9+6rUiN7ndunLvuQyKyTURODvlssXvkvMWdfq+IpIhIOvA0MCGkJXZC75ZIEblURDa73+mXRGROr2V/S0Q2iEidiPxJRFIi89cwPrAap6BZ6A6fBbwIbO81bhdwgYhsdb8fu0XkS6ELcvfp+93v3fW9viPJIvIzdz990P0/kzoC2xd1rMAZXkuADcBo4CHgEeBUnIr8M8BvRCQjZP6rgf8A8oB1wIP9LPfXQDYwDTgb+Bxwnaq2uev4TMi8VwEvqGqlOzwOSAEKgP8H3OXOfwrwAeD/icg0d95/Bj7qrmMCcAi4o1csS4FZwDL3s3NU9RmcI/w/uS1aJw30SzIx43Sc3FnZz/RbgdNwds4nAYuB74dMH4eTtwXAF4A7RGSUO+33wJdUNROYB/x9CHFd4a4nD2gD3gDecYcfBX7ea/6rgQuAE4CZwPdVtQm4CNgX0hK7L/RDIjITeBj4OpAPPAX8VUSSesVyITAVWABcO4TtMHFEVduBt3CKGNyfrwKv9Rr3ClABXAJkAdcBvxCRRXD4gPIbwPk4/1vO7rWq/8bJ84Xu9J59f/xRVXsN4QUUA41Abcjrizg7th0h880HFBgbMq4aWOi+/wPwSMi0DKALmOQOK05yBnF24nND5v0S8JL7fglQCgTc4TXAFe77c4AWIOgOZ7rLXRKyrLXAR933W4FlIdPGAx1AAlDofnZiyPS3gSvd97cBD3j997HXsOb61cCBAabvAi4OGb4AKHbf9+ReQsj0CuA09/1eN4+zei3zWuC1XuMUmO6+/wNwV8i0rwJbQ4bnA7Uhw8XAjSHDFwO7QmIs67Wuw3kM/CuwImRaACgHzglZ9mdCpv8EuNPrv5u9ovfl5tdK9/16YAZOgRw67po+Pvdn4Gvu+3uAH4VMm877/y8EaAJOCJl+OrDH62334mUtOMfmo6qaE/K6yx1/MGSeFgBV7T0utAWntOeNqjYCNTgtJ6HygCSgJGRcCU5Vjqq+hZPQZ4vIbJwk/0vIvNWq2hUaUx9x9sQ0BVjpNsfX4hQ8XTh9L3ocCHnf3Gt7jL9UA3kiktDP9AkcnZeh+Vutqp0hw6H58gmcYqPEPT17+hDi6p2/A33HIOR71keMAzli+1S1211WQcg89n0wQ/EKsNRtycxX1R3AP4Az3HHzgFdE5CIRedPtRlCL813Jc5cxgSNzOvR9PpAGrA3Zjz/jjo87VuB4a1LPG/fUVS6wr9c8VTitKKH9cybjHEn2uA/ntNNngUdVtfUY4ykFLupVvKWoavmgn3SOIIy/vAG04py27Ms+js7L3vnbJ1VdraqXAWNwjk5XuJOacHbQwOFOl8drUsj70BgHy9kjtk9ExF1WON8HY/ryBs5p2xuA1wFUtR4n125wf+4DHgN+hnMGIAfn9Ki4y9gPTAxZZmh+V+EU+SeG7MOzVTUuC28rcLx1sYgsdc/p/wfwlqqGVuO4rS8rgP8UkUxxOiJ/Awi9JPt+4GM4Rc4fjyOeO931TAEQp6NyuFfLHAQKRcRyyidUtQ7n3P0dbuf0NBFJdI8uf4LTP+X7bp7kufMOeqsAEUkSkatFJFtVO4B6nJZCcJroTxSRhW6H3duGYVO+IiITxekwfwvwJ3f8QWC025m6LyuAD4vIMhFJBL6Jc7r4H8MQk4lDqtqC043gGzj9b3q85o57BafFPhmoBDpF5CLgQyHzrgCuE5E5IpJGSP8at5XxLpw+O2MARKRARC6I3FZFL/tndGz+KkfeB6e/TpiDeQj4Ac6pqVNw+jz05as4R7a7cb4ID+GchwVAVctwOlkqR35phuqXOKe3nhORBuBNnD4+4fg/92e1iLxzHDGYKKKqP8fZ8X4fZ4dbCtyE0+ryQ5yd9QZgI04O/jDMRX8WKBaReuBG3I7yqvoezpVbLwA7cPL9eD0EPIfz/dndE6OqbsMp0na7zflHnLpS1e1uXL/GOTL+CM4tItqHISYTv17GabkMze1X3XGvqGoDzgUfK3Au9Pg0Id0OVPVp4Fc4V2DtxGkVAqf4BviuO/5N9/v1As6FIXFH3E5IZoSJyB9wOjh+f7B5w1zePThXhAzL8ozxAxEpBq5X1Re8jsWYSHBvXbAJSO7V5y3uWQuOD4hIIfBxnEtvjTHG+JiIfMw91TsK57Lwv1pxczQrcGKciPwHTvX+U1Xd43U8xhhjIu5LOKeMd+H0X/uyt+FEJztFZYwxxhjfsRYcY4wxxvhOfzfwikp5eXlaWFjodRgmRqxdu7ZKVWP6BleW82aoYj3vLefNUPWX854WOO4zNX6J8ziCu1X1xwPNX1hYyJo1a0YkNhP7RKRk8LlGluW8ibRoy3vLeRNp/eW8Z6eoxHk8/B04D7ybC1wlInO9iseYSLOcN/HGct54ycsWnMXATlXdDSAijwCXAVv6+0B1dTXr1q1j4cKFdHV1cf/997No0SIWLFhAR0cHDz74IEVFRcybN4/W1lYeeeQRlixZwpw5c2hubmbFihWcfvrpzJo1i8bGRh599FGWLl3K9OnTqaurY+XKlZx11llMmzaNQ4cO8cQTT3DOOedQWFhIVVUVTz75JMuWLWPSpElUVFTw1FNP8cEPfpCCggIOHDjAM888w4UXXsjYsWPZVVLK31e9wClnnkdK1ij2lZezdc3rzDj1bBLTsqk5WEbplrVMWvgBElIzqasoo2rHesYt+ADBlAwaK0s5tHsjYxacTSAplaaKvdSXbCF/wdkEklJoriihYe828k46l0BCEk0H99BY9h55C5cRCCbQtH83Tft2kH/yB5FAgKZ9O2nav4v8Rc4NMZvKd9BSUUzeyR8EoLFsO61VZeQtXOYMl26l7dABRi84F4CGvZtpr6ti9HznwbUNxZvoaKwhd57zENz6PRvobK4n98SlzvDudXS1NjNq7hkA1O16h+6ONkbNdh45VLdjLdrdSc4s5z6Cte+tBiBn5qnO8Pa3kEAC2TNOAeDQtjcIJCaTfcIiZ3jLPwimpJE1bSGj0pKQkrcZPXo0Z5/txPfoo48ee2ZGjq9zfndPzi89j5TMXPaVl7FlzevMPPVsEtOzqT5QTumWNX3mfEJqBg0V7+d8MDmVxoPDk/NjTrkAVfVVzqclJZBSvuaonB83bjiebDGsLOf7yPnxJzn7ecv58HM+ISBkHXx3SDnvZYFTwJEPCSujj7vmisgNOM/ooKCgoPdkT7R2dLHjYANVjW089FYJ+zsrOVRVQU5NLfctf5O9rUmM0kaKEg/x251rqNVUxgQaWZRQz92PbaReUxgbaODkhAbufXIrjZrM+EA9JyU08odnttOkSRQE6pif0MQfntlOC4lMCtRyYkIT9z67nTYSmRI4xJyEJn7/zDY6SGBqsIZZwSbufnorXQSZFqxmZrCJ5U9tQQkwPVjF9GAT//u3rQDMDFYyNdjEne7w7GAFk4KNh4fnBg8yPtjAKnd4XkIF+YFGXtzrDM9PqCBXmnm5xBk+KaGSbGnllWJn+OSEKtKlndf2OMOnJFSTLJ38Y5czfGpCNUHp5s2dzvDixBoA3t7hDJ+WeIguDbD6PWf4jMRa2jSBtduc4aWJtTRpM+9u3crMsRncGPpEpOgVsznf3tnN7spGqhvbefitvRzoqqKm0sn5P971Jntbkshxc/6OaMv5J53/pX7K+bFZyXx3dqT+2sMqZnO+s6ub4uomqpva+dPqUg68UUV1ZQXZ1XXcf7eT89ndx5jzT1vOw9ByPiUxwH+dPLS/oWeXiYvI5cAFqnq9O/xZYLGqfrW/zxQVFakX52ZVlfVldTy3+QBv76lhQ1kd7V3dAAQExmalMD47hTGZKeRmJDEqLZGslETSkxNITw6SmphASmKApIQAyQkBEgIBEoMBEoNCICAERQgGnPcBAcH5Sch7ETn8pDUAcQeOGBs6wwAkzPliRUCEjOSja3URWauqRR6E1KdYy/kt++t5dtMB3tpTw7rSWto6nZwXgTGZyYzPTmVsVjK56clOzqc6OZ+RHCQ1MUhyYpDkXjmfEHRyPTTnxV1mUCRiOR/6eT8QIDMlse9pUZT3sZTzADsONvD0pgO8taead/fW0tzedXiak/MpjMlKYXR6EqPSk8hKSSQjJYH0JCfnU5KCJAedfX1PvicGA0fkPEAwIEPKeQjJ+zjNeYCsIea8ly04ZRz5FNSJhPkk4pFyqKmde/9RzMp3yyitaSEhIMwryObaMwtZMDGbmWMzmTI6jeSEoNehmtgQ9Tnf0NrBH98o4bF3ythd2URA4MQJ2Vy9ZAoLJ+cwY0wGU/PSSUm0nDdhifqcb27v5KG39vLo2jK2HWhABGaPy+LyUyZy8uRRTB+TwQn5GaQmWc7HGi8LnNXADBGZCpQDV+I8VMxz9a0d3PXKbu59vZim9k6WTs/jn8+bwYdOHEd2at8VpDFhiNqcb27v5J7X9nDXq3uoa+lgydRcrl86jQvnjSM3Pcnr8Ezsitqcb+3o4oE3S7jz5V1UNbazaHIO/3bpiVw0fxxjMlO8Ds8MA88KHFXtFJGbgGdxLh+8R1U3exVPj3f2HuKrD71LeW0LF88fx9eWzWTWuEyvwzI+EK05v2VfPTc9/A67K5tYNnsMXz9/JvMnZnsdlvGBaM35XZWNfPWhd9myv54zp4/mzvNnUlSY63VYZph5eh8cVX0KeMrLGHqoKv/7ym5++ux2xmen8Pg/ncGiyaO8Dsv4TDTlPMD9b5bwH09uYVRaIg9dv4Qzpud5HZLxmWjL+ZXvlnHryk0kJwS4+3NFnD93rNchmQiJqTsZR9JPnt3O717axYfnj+dHn5jfb2cmY/zizpd38eOnt3HurHx+dvlJjM5I9jokYyLq4bf3cvPjG1k8NZdfXrmQ8dmpXodkIsgKHOB/X97F717axdVLJvPDj85D/Nb13JheHn57Lz9+ehuXnjSB2z+1kEDAct7425Mb9nHLyo2cOyuf5Z8rIjFoj2L0u7j/C69YU8qPnt7GR06awL9fZsWN8b9nNu0/vKP/nytOsuLG+N5rO6r4lz+to2jKKH579SlW3MSJuP4rl1Q38a9/3sTS6Xn8z+UnHb5HgTF+VVHfyrcf3cDCSTm2ozdxoa65g39ZsY5peRncfc2pdrl3HInbU1Sqyi0rN5IUDPCzy08iKcF29Mb/bvvrZto6u/n5FQttR2/iwo+f2UpNUzv3Xnuq3eYjzsTtf/XH3inn9Z3VfOei2YzLtnseGP97fstBntp4gK8tm8HUvHSvwzEm4t7cXc3Db5fyhaVTmVdgtz6IN3FZ4FQ1tvHDv22haMoorl482etwjIm4htYO/vXPm5g9LpMbzprmdTjGRFxrRxe3PL6RSbmp/Mv5M70Ox3ggLk9R3fXqbupbOvjRx+dbB0sTF+5/s4QD9a389jOLrN+NiQuPvVPG7qom/nCd9buJV3G3p2tu7+SRt0u5cN44Zoy1OxQb/+vo6ub+N0o4c/pou3mliQuqyr2vFzOvIIuzZ+Z7HY7xSNwVOI+/U05dSwefP3Oq16EYMyKe2XSA/XWtXHeG5byJD6/uqGJnRSPXnTHVbv0Rx+KqwOnuVu59fQ/zC7I5ZYodyZr4cO/re5gyOo3zZo/xOhRjRsS9r+8hLyOZS04a73UoxkNxVeC8urOKXZVNfH5poVX1Ji68u/cQ7+yt5dozCq2/mYkLuyobeXF7JZ85bTLJCdb3Jp7FVYFzz2t7yM9M5sPzJ3gdijEj4t7Xi8lMTuDyokleh2LMiPjD68UkBQNcvWSK16EYj8VNgVPR0MrL71Vy1eLJdlM/Exca2zp5ZtMBPr6ogIzkuLxg0sSZjq5unlhXzsXzx5GfaQ+PjXdx85/+pW2VAFx44jiPIzFmZLy2o4r2rm4umGc5b+LDmuJD1Ld2cqHlvCGOCpwXth5kQnYKc8bbpeEmPqzaepDMlAROLcz1OhRjRsSqrQdJCgb4wAy7NNzESYHT2tHFazurOG/OGOtcbOJCd7fy4vYKzpk1xm7sZ+LG37dVcNoJo0m3U7KGOClw3tpTQ3N7F8vmjPU6FGNGxIbyOqoa21lml4abOLG7spHdVU2W8+awuChwVm09SGpikNOnjfY6FGNGxKqtBwkInDPLmupNfPj7tgoAu9+TOcz3BY6qsmprBUtn5JGSaPdEMPFh1dYKiqbkkpOW5HUoxoyIVVsrmDU2k0m5aV6HYqKE7wuc7QcbKK9tsWZLEzf21bawZX89y+ZYzpv4UNfSweriGst5cwTfFzirtlqzpYkvPU31trM38eKV9yrp7FbLeXME3xc4q4trmDU2kzFZKV6HYuKYiPxURLaJyAYRWSkiOZFa1+riGsZlpXBCfkakVmFMVFlTXEN6UpCFk+wZg+Z9vi5wVJUNZXUsmJjtdSjGPA/MU9UFwHvAzZFaUU/O2y0RTLxYX1bHvIJsgva8NRPC1wVOeW0LNU3tVuAYz6nqc6ra6Q6+CUyMxHrqWjrYU9VkOW/iRkdXN1v211vOm6P4usDZWFYHwIKJOd4GYsyRPg883d9EEblBRNaIyJrKysohLXhzueW8iS/bDzTQ3tltOW+O4uvbPa4vqyMxKMy2xzOYESAiLwB9PQTnVlV9wp3nVqATeLC/5ajqcmA5QFFRkQ4lhvVuUT+/wI5mTXzYeLiot5w3R/J1gbOxvJbZ47JITrD735jIU9XzB5ouItcAlwDLVHVIhUu4NpbXMjk3jVHpdv8bEx82lNWRnZrIZLv/jenFt6eoejoYz7eq3kQBEbkQ+C5wqao2R2o9lvMm3mwoq7VO9aZPvi1wiqubaWjtZIE11Zvo8BsgE3heRNaJyJ3DvYLqxjbKDrVYzpu40drRxfYDDXZK1vTJt6eoNpTVAtbZ0kQHVZ0e6XVstA7GJs5s3V9PZ7da/xvTJ9+24GwsqyM5IcCMsXazMxMfNpbVIQLzCrK8DsWYEWFFvRmIbwucDWV1nDghi8SgbzfRmCOsL6tjWl46mSmJXodizIhYX1pHXkYS47PtTvXmaL7879/VrWzaV2dVvYkrG8trLedNXNlYXsv8AutgbPrmywJnT1Ujze1d1vHMxI2KhlYO1rdZzpu40dLexc6KRuZbUW/64csCZ1dlEwAzx9oN/kx82G05b+JMcXUT3QozrZ+l6YcvC5ySamdnP3m03fjJxIeenJ9iOW/iRE/OF45O9zgSE608KXBE5Kcisk1ENojIShHJGc7lF1c3k5ueRHaqdbY08aG4upnEoDAhJ9XrUIwBQEQuF5HNItItIkXDvfziaud+mXYga/rjVQvO88A8VV0AvAfcPJwLL6lusiNZE1dKqpuYlJtGMGCdLU3U2AR8HHglEgsvqW5idHoSWXbVoOmHJwWOqj6nqp3u4JvAxOFcfnFVM1PsuSQmjljOm2ijqltVdXukll9c1WytN2ZA0dAH5/PA0/1NFJEbRGSNiKyprKwcdGFtnV3sq2thip2XNXFCVd1WS8t5E5uGup8HpwXH+t+YgUTsUQ0i8gIwro9Jt6rqE+48twKdwIP9LUdVlwPLAYqKigZ9AnNpTQuqUJhnlb2JD1WN7TS1d1FoR7NmhIWznw/HUPfzrR1d7K9vta4IZkARK3BU9fyBpovINcAlwDJVHTShw7W3pudqEqvsTXw4nPN5lvNmZA22n4+UskPNzoGs7efNADx52KaIXAh8FzhbVZuHc9nFVc7iLPFNvLCcN/GmJ+etBccMxKs+OL8BMoHnRWSdiNw5XAsuqW4iMyWBUWnWs97Eh5LqJoIBocAuETdRREQ+JiJlwOnA30Tk2eFadrHdA8eEwZMWHFWdHqllF1c3Uzg63Z5NYuJGcXUzBTmpJCVEwzUDxjhUdSWwMhLLLqluJislgRw7kDUD8N0e0e6BY+KN5byJN8XVTRTm2YGsGZivCpyOrm7KDrVYs6WJKz2tlsbEi5LqZruQxAzKVwXOvtoWOrvVjmZN3KhtbqeupcNy3sSN9s5uyg41220RzKB8VeD0PJvEKnsTLyznTbwpr22hW2Gy3bnbDMJXBc77T5e1xDfxwXLexJvDV1DZfZ/MIHxV4BRXNZOaGCQ/M9nrUIwZEcVVzYjAJDuaNXGipKrnZq6W82Zgvipw9tY4V5NYz3oTrUTkWyKiIpI3HMsrqWlifFYKKYnB4VicMVGvpKaZtKQg+Rl2IGsG5qsCp7Smxc7LmqglIpOADwJ7h2uZZTUt1npj4krPft4OZM1gfFXgVDS0MjYrxeswjOnPL4DvAMP27DXLeRNvKhtaGWM5b8LgmwKnvbObQ80d1v/GRCURuRQoV9X1Ycx7g4isEZE1lZWVA85b0dBmOW/iSkVDm52eMmHx5FENkVDZ2AbAGNvZG4+IyAvAuD4m3QrcAnwonOWo6nJgOUBRUVG/rT2NbZ00t3dZzpu40d2tVDa0MSbLct4Mzj8FToNT4NjRrPGKqp7f13gRmQ9MBda7/QYmAu+IyGJVPXCs67OcN/GmtqWDzm61FhwTFt8UOBX1rQCMybRzsya6qOpGYEzPsIgUA0WqWnU8y7WcN/GmosHNeWvBMWHwTR+cnlNUdjRr4oXlvIk3Pa2WVtSbcPioBacNEcjLSPI6FGMGpKqFw7Gcinrrd2biS0/OW1FvwuGrFpzR6UkkBH2zScYMqLKxjcSgkJOW6HUoxowIu5jEDIVvqoGK+jbyrOOZiSMV9c7lsnbDMxMvKurbSEsKkp7sm5MPJoJ8U+DYzZ9MvKloaCXfct7EkYqGVmu9MWHzUYFjN38y8cVy3sSbyoY262BswuaLAkdVqWy0mz+Z+GI3PDPxptLu3G2GwBcFTm1zBx1ddvMnEz86urqpaW63nDdxxQocMxS+KHAqeu6NYEezJk5UN7ajajlv4kdLexcNbZ1W4Jiw+aTAsTu6mvhiOW/izfs5bwWOCY8vChx7Jo+JN5bzJt4cvouxXTlowuSLAufwKSrb2Zs4YTlv4k1Pzlu/MxMuXxQ4lQ128ycTX3qOZu3mliZeVFpfSzNEvihwKhra7EjWxJWKhlZGpSWSlOCLr7Axg6poaCUYEHLT7HmDJjy+2DtWNrRaXwQTV+xyWRNvKhvayMtIIhCwR5OY8PiiwKmwu1uaOGM5b+KN5bwZKl8UOJX1djRr4kuF5byJM5bzZqhivsCxmz+ZeHP40SSW8yaOWM6boYr5AqfSLpc1caa+pZP2zm4r6k3c6OpWqhutBccMTewXOI3O3S0t8U28sJw38aa6qY1utQNZMzQxX+BU1Pe04FjnMxMfLOdNvOnJ+XzLeTMEMV/gVDbaLetNfLGcN/HGct4ci5i/9e9F88YzPT+D0el28ycTHz4wI5+HvriEiaNSvQ7FmBGxaPIoHrnhNGaNy/Q6FBNDYr7Ayc9MtqrexJXc9CTOOCHP6zCMGTHZqYmcNm2012GYGBPzp6iMMcYYY3qzAscYY4wxviOq6nUMYRORSqCkj0l5QNUIhxNJftse8Gabpqhq/givc1jFUc6D/7bJq+2J6by3nI9pUZXzMVXg9EdE1qhqkddxDBe/bQ/4c5u85Mffp9+2yW/b4zU//j79tk3Rtj12isoYY4wxvmMFjjHGGGN8xy8FznKvAxhmftse8Oc2ecmPv0+/bZPftsdrfvx9+m2bomp7fNEHxxhjjDEmlF9acIwxxhhjDrMCxxhjjDG+E/MFjohcKCLbRWSniHzP63iGSkQmiciLIrJVRDaLyNfc8bki8ryI7HB/jvI61qEQkaCIvCsiT7rDMb090cRyPjpZzkeO5Xx0ivacj+kCR0SCwB3ARcBc4CoRmettVEPWCXxTVecApwFfcbfhe8AqVZ0BrHKHY8nXgK0hw7G+PVHBcj6qWc5HgOV8VIvqnI/pAgdYDOxU1d2q2g48AlzmcUxDoqr7VfUd930DTrIU4GzHfe5s9wEf9STAYyAiE4EPA3eHjI7Z7YkylvNRyHI+oizno1As5HysFzgFQGnIcJk7LiaJSCFwMvAWMFZV94Pz5QDGeBjaUN0OfAfoDhkXy9sTTSzno9PtWM5HiuV8dLqdKM/5WC9wpI9xMXndu4hkAI8BX1fVeq/jOVYicglQoaprvY7Fpyzno4zlfMRZzkeZWMn5BK8DOE5lwKSQ4YnAPo9iOWYikoiT9A+q6uPu6IMiMl5V94vIeKDCuwiH5EzgUhG5GEgBskTkAWJ3e6KN5Xz0sZyPLMv56BMTOR/rLTirgRkiMlVEkoArgb94HNOQiIgAvwe2qurPQyb9BbjGfX8N8MRIx3YsVPVmVZ2oqoU4f4+/q+pniNHtiUKW81HGcj7iLOejTKzkfEy34Khqp4jcBDwLBIF7VHWzx2EN1ZnAZ4GNIrLOHXcL8GNghYh8AdgLXO5NeMPGb9vjCcv5mOK37fGE5XxMiartsUc1GGOMMcZ3Yv0UlTHGGGPMUazAMcYYY4zvWIFjjDHGGN+xAscYY4wxvmMFjjHGGGN8xwocY4wxxviOFTjGGGOM8Z3/DxIC5qs66ZpDAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 576x288 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotT = 50\n",
    "tfp = (irf['C'] + ss['C'])/(irf['L'] + ss['L'])/((ss['C']/ss['L']))\n",
    "fig, axes = plt.subplots(2, 3, figsize=(8, 4))\n",
    "ax = axes.flatten()\n",
    "\n",
    "ax[0].plot(100*((irf['N'][0:plotT] + ss['N'])/ss['N']-1))\n",
    "ax[0].axhline(0, color='gray', linestyle=':')\n",
    "ax[0].set_title('Mass of establishments')\n",
    "ax[0].set_ylabel('pct deviation from ss')\n",
    "\n",
    "ax[1].plot(100*((irf['mu'][0:plotT] + ss['mu'])/ss['mu']-1))\n",
    "ax[1].axhline(0, color='gray', linestyle=':')\n",
    "ax[1].set_title('Markup')\n",
    "\n",
    "ax[2].plot(100*((tfp[0:plotT] - 1)))\n",
    "ax[2].axhline(0, color='gray', linestyle=':')\n",
    "ax[2].set_title('TFP')\n",
    "\n",
    "ax[3].plot(100*((irf['L'][0:plotT] + ss['L'])/ss['L']-1))\n",
    "ax[3].axhline(0, color='gray', linestyle=':')\n",
    "ax[3].set_title('Employment')\n",
    "\n",
    "ax[4].plot(100*((irf['C'][0:plotT] + ss['C'])/ss['C']-1))\n",
    "ax[4].axhline(0, color='gray', linestyle=':')\n",
    "ax[4].set_title('Consumption')\n",
    "\n",
    "ax[5].plot(100*((irf['w'][0:plotT] + ss['w'])/ss['w']-1))\n",
    "ax[5].axhline(0, color='gray', linestyle=':')\n",
    "ax[5].set_title('Wage')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "60d1b908",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = [100*((irf['N'][0:plotT] + ss['N'])/ss['N']-1),\n",
    "       100*((irf['mu'][0:plotT] + ss['mu'])/ss['mu']-1),\n",
    "       100*((tfp[0:plotT] - 1)),\n",
    "       100*((irf['L'][0:plotT] + ss['L'])/ss['L']-1),\n",
    "       100*((irf['C'][0:plotT] + ss['C'])/ss['C']-1),\n",
    "       100*((irf['w'][0:plotT] + ss['w'])/ss['w']-1)]\n",
    "data = np.array(data)\n",
    "np.savetxt(\"ces_irf.csv\", np.transpose(data), delimiter=\",\", fmt = \"%g\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
